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ABSTRACT 

By  providing  highly  realistic  simulations  of  sound 
propagation  through  complex  atmospheric  and  terrain  en¬ 
vironments,  finite-difference  time-domain  (FDTD)  tech¬ 
niques  can  potentially  reduce  development  time  and  im¬ 
prove  the  battlefield  performance  of  acoustic  sensors.  In 
this  paper,  we  summarize  recent  progress  in  improving 
two  key  aspects  of  acoustic  FDTD  calculations  for  the  at¬ 
mosphere:  (1)  development  of  a  rigorous  implementation 
of  sound  propagation  in  a  moving,  inhomogeneous  fluid, 
and  (2)  formulation  and  numerical  implementation  of  time- 
domain  methods  for  handling  sound  interactions  with  par¬ 
tially  reflecting  ground  surfaces.  The  new  techniques  are 
illustrated  with  highly  detailed  calculations  of  sound  prop¬ 
agation  through  simulated,  dynamic  atmospheric  turbu¬ 
lence  fields  and  over  a  porous  ground  surface  with  viscous 
and  thermal  relaxation  mechanisms. 

1.  INTRODUCTION 

Acoustic  sensors  are  expected  to  play  a  key  role  in 
the  Army’s  Future  Force  by  providing  rapidly  deployable, 
networked  surveillance  over  wide  areas.  It  is  well  known 
that  the  performance  of  acoustic  sensors  is  affected  by 
complex  sound  propagation  phenomena  occurring  in  out¬ 
door  settings,  such  as  reflections  from  trees  and  buildings, 
ground  interactions,  scattering  by  turbulence,  refraction 
by  wind  and  temperature  gradients,  and  diffraction  over 
hills.  The  expense  and  difficulty  of  performing  comprehen¬ 
sive,  controlled  field  experiments  outdoors,  combined  with 
the  rapid  schedule  for  fielding  the  Future  Force,  makes 
necessary  novel  approaches  to  advanced  sensor  develop¬ 


ment.  A  realistic  simulation  of  complex  environmental  ef¬ 
fects  on  propagating  sound,  when  combined  with  detailed 
source  signature  and  high-resolution  atmospheric  and  ter¬ 
rain  inputs,  could  be  immensely  valuable  in  shortening 
the  requirements-to-deployment  cycle  for  acoustic  sensors, 
developing  robust  signal  information  processing  stategies, 
and  improving  sensor  utilization  doctrine. 

Most  current  numerical  methods  for  outdoor  sound 
propagation,  such  as  the  fast  field  program  and  parabolic 
equation  (Salomons,  2001),  are  incapable  of,  or  poorly 
suited  to,  simulating  all  of  the  propagation  phenomena 
mentioned  in  the  preceeding  paragraph.  Complex,  moving 
source  distributions,  such  as  maneuvering  ground  vehicles, 
are  also  difficult  to  incorporate.  These  shortcomings  can 
potentially  be  overcome  with  finite-difference,  time-domain 
(FDTD)  techniques,  which  have  become  popular  for  elec¬ 
tromagnetic  and  seismic  wave  propagation.  But  there 
are  significant  drawbacks  to  FDTD  techniques  that  have 
so  far  prevented  their  widespread  use  for  outdoor  sound 
propagation.  Foremost  among  these  is  that  they  are  very 
computationally  intensive  when  applied  to  the  frequency 
range  (a  few  hundred  Hz  and  lower)  and  spatial  scales  (a 
few  km  or  less)  of  Army  tactical  applications  for  acoustic 
sensors.  Fortunately,  the  current  generation  of  parallel¬ 
processing  computers  now  makes  FDTD  calculations  vi¬ 
able.  Some  other  difficulties,  particular  with  the  sound 
propagation  problem,  are  (1)  incorporating  the  dynami¬ 
cally  moving  (windy  and  turbulent)  atmospheric  propaga¬ 
tion  medium,  and  (2)  formulating  time-domain  techniques 
for  the  acoustic  interactions  with  the  ground,  including  ab¬ 
sorption  and  dispersion  characteristics  of  porous  materials 
such  as  soils.  Successfully  addressing  these  two  problems  is 
key  to  making  acoustic  FDTD  simulation  useful  for  Army 
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applications.  These  problems  are  the  focus  of  this  paper. 

Incorporation  of  a  dynamically  moving  atmosphere  is 
discussed  in  Section  2.  There  the  time-domain  differential 
equations  needed  to  propagate  the  sound  fields  through 
a  propagation  medium  in  motion  are  described  and  some 
of  the  numerical  issues  involved  in  their  implementation 
are  addressed.  An  example  is  provided  with  inputs  to  an 
FDTD  calculation  provided  by  an  atmospheric  large-eddy 
simulation  (LES).  Section  3  describes  acoustic  interactions 
with  porous  materials  such  as  the  ground  and  their  im¬ 
plementation  in  time-domain  calculations.  Due  to  the  his¬ 
torical  emphasis  in  acoustics  on  frequency-domain  compu¬ 
tational  methods,  our  time-domain  analysis  represents  a 
fundamentally  new  approach. 

2.  ACOUSTIC  FDTD  CALCULATIONS  IN  A 
MOVING  MEDIUM 

For  most  terrestrial  problems  involving  electromagnet¬ 
ics  or  seismics,  the  propagation  medium  (atmosphere  or 
Earth)  does  not  move  or  otherwise  change  significantly  as 
waves  propagate  between  sources  and  receivers  of  inter¬ 
est.1  The  same  cannot  be  said,  however,  for  sound  propa¬ 
gation  through  the  atmosphere.  Propagation  times  and  at¬ 
mospheric  variations  both  occur  in  seconds  or  less  for  most 
scenarios  of  tactical  interest  to  the  Army.  Wind  Mach  num¬ 
bers  (the  ratio  of  the  wind  speed  to  the  sound  speed)  are 
commonly  as  high  as  1/50  in  the  near-Earth  atmosphere 
and  1/5  in  a  stratospheric  jet  stream.  To  be  useful  in  Army 
applications,  numerical  methods  for  sound  propagation  in 
the  atmosphere  must  account  for  the  effects  of  wind,  tur¬ 
bulence,  and  other  atmospheric  disturbances. 

In  this  section,  we  consider  various  aspects  of  including 
a  moving  propagation  medium  in  acoustic  FDTD  calcula¬ 
tions.  An  accurate  treatment  has  required  derivation  of  a 
new  set  of  coupled  state  equations,  together  with  numerical 
methods  for  solving  it. 

2.1  Coupled  Equation  Set 

The  wave  equation,  which  is  the  starting  point  for 
most  sound  propagation  calculations,  is  a  second-order  par¬ 
tial  differential  equation  in  both  time  and  space.  FDTD 
techniques,  however,  are  most  readily  applied  to  first-order 
partial  differential  equations  (that  is,  a  state  equation  set). 
Furthermore,  most  solutions  of  the  wave  equation  have 
been  based  on  one-way  approximations,  in  which  the  en¬ 
ergy  is  propagated  in  only  one  direction,  and  on  effective 
sound-speed  approximations,  in  which  the  sound  speed  is 
taken  to  be  the  actual  sound  speed  plus  the  component  of 
the  wind  velocity  in  the  direction  of  propagation.  Wave 
equations  in  a  moving  medium  that  do  not  use  these  ap¬ 
proximations  are  considerably  more  complicated  than  the 
wave  equation  for  a  stationary  medium  (Ostashev,  1997). 
Fortunately,  the  switch  to  first-order  equations  facilitates 
correct  handling  of  the  wind  velocity  field.  The  following 
coupled,  first-order  equations  for  the  acoustic  pressure  p 


1  There  are  some  notable  exceptions  to  this  statement,  such 
as  Doppler  shifts  in  signals  from  clear-air  radars. 


and  acoustic  particle  velocity  w  involve  no  one-way  or  ef¬ 
fective  sound-speed  approximations,  and  therefore  provide 
an  appropriate  starting  point  for  accurate  FDTD  calcula¬ 
tions  in  a  moving  atmosphere  (Ostashev  et  ah,  2004): 

^  =  -  (v  •  V)  p  -  pc2  V  •  w +pc2Q,  (1) 

^  = -(w- V)v- (v- V)w- ^  +  T  (2) 

Here,  p  is  medium  density,  c  is  the  adiabatic  speed  of  sound, 
and  v  is  the  wind  velocity.  The  quantities  F  and  Q  repre¬ 
sent  sources:  the  former  is  a  force  acting  on  the  medium, 
whereas  the  latter  is  a  mass  source.  Bold  symbols  represent 
vectors.  The  terms  involving  (v  •  V),  which  are  particular 
to  the  moving  medium,  are  called  the  advective  terms.  Nu¬ 
merical  issues  aside,  the  source  terms  and  atmospheric  field 
variables  can  be  arbitrary  functions  of  space  and  time. 

Equations  (1)  and  (2)  were  derived  in  Ostashev  et  al. 
(2004)  from  the  full  fluid  dynamical  equations  using  the 
conventional  linear  acoustics  approximation,  namely  that 
the  sound  wave  is  a  small  perturbation  to  the  background 
state  of  the  medium.  Additionally,  the  sound  waves  are 
assumed  to  be  uninfluenced  by  the  following:  (1)  diver¬ 
gence  in  the  atmospheric  flow,  (2)  the  background  pres¬ 
sure  gradient,  and  (3)  the  Coriolis  force.  The  first  of  these 
implies  that  turbulent  fluctuations  in  the  flow  have  a  low 
Mach  number,  which  is  very  reasonable  for  the  problems 
considered  here.  The  second  may  be  considered  a  distin¬ 
guishing  property  between  sound  and  gravity  waves.  The 
third  is  valid  for  even  for  very  low  infrasonic  frequencies. 
An  FDTD  code  based  on  Eqs.  (1)  and  (2)  is  both  more 
general  and  accurate  than  most  current  sound  propagation 
formulations,  despite  its  comparative  simplicity. 

2.2  Computational  Considerations 

Typically,  finite-difference  solutions  for  wave  propaga¬ 
tion  in  a  nonmoving  medium  use  a  grid  that  is  staggered  in 
space  and  time  (Botteldooren,  1994;  Graves,  1996).  Each 
acoustic  particle  velocity  component  is  explicitly  calculated 
on  spatial  grid  nodes  shifted  by  one-half  of  the  internode 
spacing,  relative  to  the  explicit  acoustic  pressure  nodes,  in 
the  direction  of  the  velocity  component.  The  particle  veloc¬ 
ities  and  pressures  are  advanced  on  alternating  time  steps. 
The  finite-difference  stencil  corresponding  to  the  pressure 
advancement  in  this  approach  is  illustrated  in  Figure  1.  It 
happens,  however,  that  this  staggered  “leap-frog”  method¬ 
ology  of  alternately  marching  the  solution  in  time  cannot 
be  applied  directly  to  Eqs.  (1)  and  (2).  Evaluation  of  the 
advective  terms  on  the  right-hand  sides  of  these  equations 
requires  knowledge  of  the  pressure  and  particle  velocity 
fields  at  time  steps  where  they  are  not  explicitly  available. 
Therefore,  we  have  developed  alternative  approaches  based 
on  unstaggered  temporal  grids  (Ostashev  et  al.,  2004; 
Wilson  and  Liu,  2004)  and  staggered  temporal  grids  span¬ 
ning  multiple  time  steps  (Symons  et  al.,  2003).  (In  either 
case,  a  staggered  spatial  grid  is  still  used.)  The  latter,  stag¬ 
gered  approach  is  shown  schematically  in  Figure  2.  This 
finite-difference  stencil  provides  fully  centered  spatial  finite 
differences  in  the  context  of  the  moving  medium  problem. 
As  discussed  in  Wilson  and  Liu  (2004),  this  approach  and 


2 


in  which  a  source  emits  frequencies  at  250  Hz  and  lower. 
For  a  sound  speed  of  340  ms-1,  the  minimum  wavelength 
is  340/250  =  1.36  m,  and  the  spatial  grid  interval  therefore 
0.17  m.  Assuming  the  dimensions  of  the  computational 
domain  are  500  m  in  each  horizontal  direction  and  50  m  in 
the  vertical,  about  2.5  x  109  grid  nodes  are  required.  Such 
intensive  computational  problems  can  be  tackled  only  on 
very  large,  parallel  processing  computers. 

2.3  High  Mach  Number  Tests  of  the  Solution 
Technique 


Fig.  1:  Standard  0(4,2)  (fourth  order  in  space,  second 
order  in  time)  staggered  leap-frog  finite-diference  stencil 
for  updating  the  acoustic  pressure  solution  in  a  non-moving 
medium.  For  simplicity,  the  stencil  is  shown  with  only  one 
spatial  dimension.  Time  is  the  vertical  axis  and  space  is 
horizontal,  with  corresponding  grid  intervals  At  and  Ax, 
respectively.  The  dashed  lines  intersect  at  the  pivot  point, 
which  is  the  location  on  which  the  finite-difference  approx¬ 
imations  are  centered. 


Fig.  2:  New  0(4,2)  staggered  finite-diference  stencil 
for  updating  the  acoustic  pressure  solution  in  a  moving 
medium.  Two  time  levels  must  be  stored  in  comparison  to 
the  single  time  level  for  the  non-moving  medium. 


several  others  can  yield  highly  accurate  results,  although 
some  efficiency  in  memory  usage  and/or  calculation  time 
is  lost  in  comparison  to  the  customary  leap-frog  solution 
for  a  nonmoving  medium. 

The  size  and  dimensions  of  the  computational  domain 
in  an  acoustic  FDTD  simulation  depend  on  the  propaga¬ 
tion  geometry  of  interest  and  the  memory  available.  Ar¬ 
tificial  sound  absorbing  layers  are  usually  placed  around 
the  sides  and  corners  of  the  domain  to  prevent  unwanted 
numerical  reflections.  The  lower  surface  is  normally  taken 
as  a  rigid  plate,  although  more  sophisticated  and  realistic 
partially  absorbing  ground  models  are  under  development 
as  described  in  Section  3.  The  spacing  between  grid  points 
in  acoustic  is  driven  by  the  smallest  acoustic  wavelength 
of  interest.  For  fourth-order  spatial  finite  differencing,  a 
typical  grid  spacing  would  be  1/8  of  the  wavelength.  To 
illustrate  the  memory  requirements,  consider  a  simulation 


We  now  consider  2D  FDTD  calculations  for  the  case 
of  a  stationary  source  in  a  uniform,  high  Mach  number 
wind.  Although  this  is  a  simple  case,  it  is  very  useful  for 
testing  the  fidelity  of  the  numerical  solution  method,  since 
an  analytical  solution  is  known.  For  this  series  of  calcula¬ 
tions,  a  mass-type  source  was  used  that  consisted  of  a  10- 
cycle,  100-Hz  signal  with  a  cosine  taper  function  applied  to 
the  three  periods  at  the  beginning  and  end.  The  tapering 
alleviates  numerical  dispersion  of  high  frequencies,  which 
can  become  evident  when  there  is  an  abrupt  change  in  the 
source  emission. 

Figure  3  illustrates  the  pressure  field  calculation  for 
a  uniform  Mach  0.3  wind.  The  field  is  shown  at  0.11  s, 
or  0.01  s  after  the  source  has  been  turned  off.  Note  that 
the  distance  between  wave  fronts  is  smaller  upwind  than 
downwind,  due  to  the  role  of  the  wind  in  determining  the 
overall  propagation  speed.  Although  not  as  clearly  evi¬ 
dent,  the  pressure  amplitude  is  higher  upwind  than  down¬ 
wind.  The  azimuthal  dependence  of  \p(r,  a,  M)/p(r,  0,  0)| 
at  a  distance  r  =  (10/7t)A  (where  A  is  the  wavelength,  a 
the  azimuthal  angle  measured  from  upwind,  and  M  is  the 
Mach  number)  is  plotted  and  compared  to  the  theoretical 
prediction  in  Figure  4  for  M  —  0,  0.3,  and  0.6.  Two  cal¬ 
culated  curves  for  each  value  of  M  are  shown,  one  for  a 
low-resolution  run  with  800x800  grid  points  and  a  time 
step  At  =  0.145  ms,  and  the  other  for  a  high-resolution 
run  with  1600x1600  grid  points  and  At  =  0.0362  ms.  For 
M  =  0.3,  both  grid  resolutions  yield  nearly  exact  agree¬ 
ment  with  theory.  Agreement  with  theory  at  M  =  0.6  is 
very  good  for  the  high-resolution  run,  although  the  low- 
resolution  run  substantially  underpredicts  the  upwind  am¬ 
plitude.  High  spatial  resolution  is  needed  at  high  Mach 
numbers  because  of  the  shortening  of  the  wavelength  in 
the  upwind  direction. 

2.4  Simulations  with  Dynamic  Atmospheric 
Models 

The  time-domain  nature  of  acoustic  FDTD  simulation 
makes  it  a  natural  approach  for  coupling  with  dynamic 
atmospheric  models  such  as  global  and  mesoscale  numeri¬ 
cal  weather  predictions  (NWP)  and  large-eddy  simulations 
(LES)  of  the  atmospheric  boundary-layer  (ABL)  and  flow 
in  urban  terrains.  In  this  section,  we  present  examples  of 
FDTD  simulations  of  sound  propagating  through  ABL  tur¬ 
bulence  fields  generated  by  LES.  The  LES  we  use  is  a  par¬ 
allelized  implementation  of  the  physical  models  and  code 
described  in  Sullivan  et  al.  (1994).  Two  stability  cases 
are  considered.  The  first  is  for  a  neutral  ABL  (Figure  5) 
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Fig.  3:  Wavefronts  of  the  sound  pressure  due  to  a  100-Hz 
point  source  located  at  x  =  0  and  y  —  0.  A  uniform  flow 
with  Mach  number  0.3  moves  from  left  to  right. 


Mach  0.3,  Time  -  0*110  s 


Fig.  4:  Comparison  of  the  angular  dependence  of  the 
pressure  amplitude  from  FDTD  calculations  to  theory.  The 
downwind  direction  is  0°  and  the  upwind  direction  is  180°. 
Calculations  for  grids  at  two  different  resolutions  and  three 
different  Mach  numbers  are  shown. 


and  the  second  for  a  buoyantly  unstable  ABL  (Figure  6). 
Although  only  the  vertical  velocity  component  is  shown  in 
these  figures,  all  three  velocity  components  and  temper¬ 
ature  were  used  as  input  to  the  acoustic  FDTD  calcula¬ 
tion.  The  sound  speed  c  in  Eqs.  (1)  and  (2)  is  known  from 
the  absolute  temperature  T  according  to  the  relationships 
c  =  20.02 \fT\  the  density  follows  from  the  ideal  gas  law. 

The  acoustic  calculation  was  performed  in  3D  with 
901  x  901  x  603  grid  nodes  spaced  at  1  m.  The  time  step 
was  0.25  ms  and  the  solution  was  advanced  over  16001  time 
steps.  The  sound  source  was  a  20-Hz,  mono- frequency, 


mass-type  source.  Execution  required  10  hours  on  a  clus¬ 
ter  built  from  100  Compaq  ES45  processors,  making  this 
run  likely  the  single  most  computationally  intensive  and 
detailed  calculation  of  sound  propagation  in  the  ABL  yet 
performed. 

The  middle  panels  in  Figures  5  and  6  show  snapshots 
of  the  calculated  sound  fields.  Distortions  to  the  propagat¬ 
ing  wavefronts  are  not  easily  discernable  in  these  images  as 
presented.  For  each  case  we  also  propagated  sound  fields 
through  horizontally  averaged  LES  fields  (i.e.,  the  mean 
vertical  profiles  only).  The  difference  between  the  sound 
fields,  with  and  without  the  horizontal  LES  variability,  is 
shown  in  the  lower  panels  of  Figures  5  and  6.  The  dis¬ 
tortions  to  wavefront  shape  and  amplitude  are  easily  dis¬ 
cernable  in  these  difference  images.  Such  distortions  cause 
fluctuations  in  apparent  bearings  of  targets  derived  from 
acoustic  sensor  arrays  and  are  the  limiting  factor  in  array 
performance  when  the  signal-to-noise  ratio  is  high.  The  ca¬ 
pability  to  realistically  simulate  turbulence  effects  enables 
virtual  testing  of  acoustic  beamforming  systems  being  de¬ 
veloped  for  the  Army’s  Future  Force. 


3.  TIME-DOMAIN  MODELING  OF  AN 
ABSORPTIVE  GROUND 


Typically,  in  frequency-domain  acoustic  calculations, 
the  ground  interaction  is  characterized  using  a  localized 
impedance  function  Z,  defined  as 


P  (z  =  0,  u 
Wn  (Z  0,  OJ) 


in  which  P  and  Wn  are  the  Fourier  transforms  of  the  pres¬ 
sure  and  particle  velocity  normal  to  the  boundary,  respec¬ 
tively,  z  —  0  is  the  position  of  the  boundary,  uo  —  2i r/,  and 
/  is  the  frequency.  In  this  formulation,  the  acoustic  pres¬ 
sure  is  analogous  to  voltage  in  an  electrical  system,  while 
the  particle  velocity  is  analogous  to  current.  The  linear 
ground  boundary  condition  between  the  pressure  and  par¬ 
ticle  velocity  generally  makes  frequency-domain  handling 
of  the  ground  interaction  a  simple  task. 

Alternatively,  one  might  explicitly  calculate  the  wave 
propagation  in  the  ground  as  well  as  in  the  air.  In  that 
case,  the  normal  methodology  is  to  use  complex  opera¬ 
tors  (frequency-dependent  values  with  real  and  imaginary 
parts)  to  represent  the  bulk  material  properties,  specifi¬ 
cally  the  density  and  bulk  modulus  within  the  ground.  The 
imaginary  parts  of  the  operators  relate  to  wave  attenuation 
and  dispersion  in  the  medium.  Implementation  of  complex¬ 
valued  material  properties  in  most  computer  programming 
languages  is  little  more  difficult  than  real-valued  material 
properties. 

The  impedance  ground  boundary  condition  and  bulk 
material  methods  each  have  their  own  relative  merits.  The 
former  is  highly  computationally  efficient:  the  wavefield  in 
the  ground  does  not  need  to  be  calculated  or  stored.  On 
the  other  hand,  the  bulk  material  method  is  called  for  when 
one  is  interested  in  the  waves  in  the  ground  or  they  must 
be  explicitly  calculated  due  to  reflections  from  non-uniform 
ground  properties,  such  as  when  a  thin  layer  of  soil  overlies 
rock.  The  bulk  material  method  is  also  easier  to  apply 
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Neutral  LES:  V  Wind 


Unstable  LES:  V  Wind 


200  400  600  800 


Fig.  5:  Top:  Vertical  wind  field  from  an  LES  for  a  neu¬ 
tral  ABL  ( Zi/L0  =  —  oo).  Only  a  partial  cross  section  of  the 
full  LES  is  shown.  The  turbulent  structures  are  generated 
primarily  by  wind  shear.  Middle:  Acoustic  FDTD  calcula¬ 
tion  of  waves  from  a  constant-frequency  source  propagating 
through  the  neutral  LES  fields.  Bottom:  Difference  field, 
multiplied  by  10,  between  sound  field  calculated  with  tur¬ 
bulence  (LES  fields)  and  without  turbulence  (horizontally 
averaged  LES  fields). 


when  the  surface  is  not  parallel  to  one  of  the  coordinate 
directions. 

Regardless  of  whether  a  ground  boundary  condition  or 
bulk  material  approach  is  appropriate,  time-domain  proce¬ 
dures  for  handling  attenuation  and  dispersion  in  a  porous 
ground  material  are  not  nearly  so  well  established  as  the 
frequency-domain  procedures.  Fundamentally,  the  issue  is 
the  reaction  time  of  the  porous  material,  which  requires  the 


Fig.  6:  Same  as  Figure  5,  except  for  a  buoyantly  unstable 
ABL  (zi/Lo  =  -6). 


history  of  the  signal  to  be  stored  over  an  interval  of  time. 
If  this  interval  is  long  compared  to  the  computational  time 
step,  a  substantial  amount  of  computer  memory  may  be 
required.  In  the  remainder  of  this  section,  we  describe  ap¬ 
proaches  to  time-domain  modeling  of  the  ground,  based 
both  on  the  bulk  material  and  impedance  ground  bound¬ 
ary  condition  methods. 

3.1  Bulk  Material  Method 

Although  dozens  of  frequency-domain  models  have 
been  developed  for  sound  propagation  in  porous  materi¬ 
als,  we  have  found  that  one  based  on  explicitly  modeling 
the  viscous  and  thermal  diffusion  in  the  pores  as  relax¬ 
ation  processes  (Wilson,  1993)  translates  most  naturally  to 
the  time  domain.  Other  approaches  that  have  been  exam- 
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ined  are  either  non-causal  (Berthelot,  2001)  or  limited  in 
their  frequency  range  of  applicability  (Zwikker  and  Kosten, 
1949;  Fellah  and  Depollier,  2000),  and  in  any  case  do  not 
have  closed-form  analytical  solutions  in  the  time  domain. 
A  new  set  of  integro-differential  equations  (Wilson  et  ah, 
2004)  to  be  discussed  in  this  section  remedies  all  of  these 
shortcomings  with  no  appreciable  loss  in  accuracy. 

The  physics  underlying  the  relaxation  model  can  be 
understood  in  a  straight-forward  manner.  When  a  pressure 
gradient  is  introduced  across  a  porous  material  (by  a  sound 
wave  or  other  mechanism),  air  is  accelerated  through  the 
pores.  The  acceleration  is  counteracted  by  viscous  drag 
forces  exerted  by  the  solid  frame  material.  The  air  flow 
gradually  relaxes  to  a  steady  state  in  which  the  pressure 
gradient  and  viscous  forces  balance.  Similarly,  if  the  air  in 
the  pores  is  heated  or  cooled,  the  temperature  in  the  pores 
gradually  relaxes  to  a  thermal  equilibrium  with  the  frame 
material  as  heat  is  conducted  between  the  solid  and  air.2 

The  new  set  of  integro-differential  equations  (Wilson 
et  ah,  2004)  describing  these  processes  in  the  time  domain 
is: 


w  ^  <9w  ^  1  p  w  ( t' )  /rv  +  <9w  (tf)  /dt' 

Tv  +  dt  +  ,/7rTV  Too  -  t 7 

x  exp  (—  — — t )  dt'  =  _p'ooVp;  (3) 


and 


o  dp  Poo  (7  - 1)  r*  dpj/y&P 

p°° dt +  Jirr;  ]_„ 


exp  (—  - — — )  dt'  =  — V  •  w.  (4) 


Here,  rv  is  the  time  constant  for  the  viscous  relaxation 
process  and  re  for  the  thermal  (entropy)  relaxation  process. 
Additional  symbols  are:  T4o  =  Q/pq2 ,  the  effective  specific 
volume  at  high  frequency,  /3oo  =  Q/P^y,  the  effective  (adi¬ 
abatic)  compressibility  at  high  frequency,  £4,  the  porosity 
(void  fraction),  q ,  the  tortuosity  (a  measure  of  the  oblique¬ 
ness  of  the  pores),  P,  the  ambient  air  pressure,  and  7,  the 
ratio  of  specific  heats  for  air. 


The  integrals  in  Eqs.  (3)  and  (4)  represent  convolu¬ 
tions  of  the  propagating  wavefield  with  the  relaxational 
response  function  s  (t)  =  1/y/irTt  exp  (—t/r)  H  (t)  [where 
H  ( t )  is  the  Heaviside  function  equal  to  0  for  t  <  0  and 
1  for  t  >  0].  When  the  relaxation  time  is  very  short  in 
comparison  to  any  changes  in  the  propagating  wavefield, 
the  response  functions  can  be  approximated  as  unit  im¬ 
pulses  and  the  following  much  simpler  set  of  equations  is 


obtained: 


2 

—  w  + 

Tv 


3  (tw 
2  dt 


—  VooVp, 


^7S="v'w- 


The  first  of  these  equations  describes  a  simple  balance 
among  viscous,  inertial,  and  pressure  gradient  forces, 
whereas  the  second  implies  that  the  propagation  is  isother¬ 
mal  (that  is,  the  acoustic  period  is  long  enough  that  the 
temperature  in  the  pores  stays  in  equilibrium  with  the 
frame  material).  With  notational  changes,  these  equations 
are  the  same  as  a  well  known  set  of  phenomenological  equa¬ 
tions  suggested  by  Zwikker  and  Kosten  (1949)  that  has 
more  recently  been  used  for  acoustic  FDTD  calculations 
(Salomons  et  ah,  2002). 

As  discussed  in  Section  2.2,  in  FDTD  simulation  the 
wavefield  variables  are  calculated  at  discrete  time  steps. 
The  wavefield  variables  should  vary  only  a  small  amount 
between  the  time  steps.  Using  this  assumption,  but  avoid¬ 
ing  any  additional  ones  pertaining  to  the  relative  values  of 
the  time  step  and  relaxation  times,  we  can  readily  evaluate 
the  integrals  in  Eqs.  (3)  and  (4)  in  a  closed  form  involving 
error  functions.  Although  the  details  are  not  given  here, 
we  have  succeeded  in  deriving  such  a  solution  and  imple¬ 
menting  it  in  a  2D  FDTD  code. 

An  example  calculation  involving  a  porous  ground  sur¬ 
face  is  shown  Figures  7-10.  The  domain  is  configured  with 
a  source  in  the  center,  a  rigid  barrier  20  m  to  the  right  of 
the  source,  a  20-m  thick  porous  ground  layer  at  the  bottom 
of  the  domain,  and  a  20-m  thick  artificial  absorbing  layer 
at  the  top.  The  porous  material  parameters  for  the  ground 
are  a  —  1000  Pa-s-m_2,  q  =  1.8,  and  Q  =  0.5,  which  are 
reasonable  for  an  “acoustically  soft”  ground  such  as  snow 
or  coarse  gravel.  The  source  in  the  simulation  emits  10 
cycles  of  a  100-Hz  sine  wave,  as  described  in  Section  2.2. 
No  wind  flow  is  present  in  this  calculation. 

In  Figure  7,  corresponding  to  0.058  s  after  source  ini¬ 
tiation,  the  sound  waves  are  just  beginning  to  impinge  on 
the  barrier.  By  the  time  of  the  second  snapshot  (Figure  8, 
taken  at  0.116  s)  the  wavetrain  is  just  starting  to  penetrate 
the  ground.  A  diffracted  wave  curls  over  the  right  side  of 
the  barrier  and  a  reflection  is  evident  to  the  right.  At  0.151 
s  (Figure  9),  the  reflection  off  the  barrier  is  fully  formed  and 
a  partial  reflection  from  the  ground  is  also  becoming  evi¬ 
dent.  The  sound  in  the  ground  has  a  shorter  wavelength 
than  in  the  air  and  is  rapidly  attenuated.  In  Figure  10, 
taken  at  0.250  s,  the  initial  wave  and  barrier  reflection  are 
distinct  and  propagating  leftward,  while  the  barrier  dif¬ 
fraction  propagates  to  the  right.  Weak  ground  reflections 
associated  with  the  initial  wave  and  barrier  reflection  are 
also  visible.  Waves  in  the  ground  refract  strongly  toward 
the  surface  normal,  as  is  consistent  with  Snell’s  law. 

3.2  Ground  Boundary  Condition  Method 

As  mentioned  earlier,  a  ground-boundary  condition 
formulation  can  provide  a  substantial  savings  in  computa¬ 
tional  memory  and  processing  time  because  the  wavefield 
need  not  be  explicitly  calculated  in  the  ground.  Based  on 
the  viscous  and  thermal  relaxation  concept,  Wilson  (1993) 
previously  derived  a  full-frequency  equation  for  the  im¬ 
pedance  of  a  porous  ground  surface: 


2  Typically,  the  heat  capacity  of  the  frame  material  is  so  large 
compared  to  the  air  that  the  frame  essentially  remains  at  a  fixed 
termperature  while  the  temperature  of  the  air  in  the  pores  at¬ 
tempts  to  attain  equilibrium. 


(7) 
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Time  =  0*058  s  Time  =  0*250  s 


Fig.  7:  Propagation  above  a  soft  porous  ground  in  the 
presence  of  a  rigid  barrier,  0.058  s  after  source  initiation. 
The  100-Hz  source  is  located  at  the  middle  of  the  domain, 
and  an  absorbing  layer  is  present  at  the  top  to  eliminate 
numerical  reflections. 
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Fig.  8:  Propagation  above  a  soft  porous  ground  in  the 
presence  of  a  rigid  barrier,  0.116  s  after  source  initiation. 


Time  =  0*116  s 


Time  =  0*151  s 


Fig.  9:  Propagation  above  a  soft  porous  ground  in  the 
presence  of  a  rigid  barrier,  0.151  s  after  source  initiation. 


In  principle,  an  inverse  Fourier  transform  could  be  ap¬ 
plied  to  this  equation  to  determine  a  time-domain  ground 
boundary  condition.  However,  a  closed-form  inverse  trans¬ 
formation  is  not  known.  Fortunately,  we  have  found  that 
the  following  much  simpler  equation,  derived  from  a  fac- 


Fig.  10:  Propagation  above  a  soft  porous  ground  in  the 
presence  of  a  rigid  barrier,  0.250  s  after  source  initiation. 


torization  of  (7),  behaves  nearly  identically: 


2M  = 


pcq  1  —  i(jJTz 

Q  V  —iwTz  5 


(8) 


where  rz  =  (y/2 )rv.  This  approximation  is  compared  to 

(7)  in  Figure  11.  Quite  interestingly,  one  can  show  that 

(8)  has  the  same  functional  form  as  an  impedance  derived 
from  the  low-frequency  approximations  (5)  and  (6).  Phys¬ 
ically,  this  suggests  that  the  convolution  terms  in  (3)  and 
(4),  which  cause  dissipation  of  propagating  acoustic  en¬ 
ergy  at  higher  frequencies,  do  not  significanly  affect  the 
impedance.  The  inverse  transform  of  Eq.  (8)  can  be  found 
in  tables.  It  is 


where  To  and  Ii  are  the  modified  Bessel  functions  of  zeroth 
and  first  order,  respectively. 

The  FDTD  calculation  for  propagation  with  a  barrier 
and  soft  porous  ground  from  Section  3.1  is  repeated  in  Fig¬ 
ure  12,  except  that  the  ground  interaction  in  this  instance 
is  calculated  with  Eq.  (9).  Although  the  propagation  in 
the  ground  is  not  explicitly  calculated,  the  appearance  of 
the  ground  reflections  is  nearly  the  same  as  in  Figure  10. 

4.  CONCLUSION 

Improvements  in  acoustic  FDTD  techniques,  and  the 
increasing  capabilities  of  modern  parallel-processing  com¬ 
puters,  make  possible  highly  detailed  calculations  of  sound 
propagation  through  the  atmosphere  with  real-world  fea¬ 
tures  such  as  complex  reflecting  surfaces,  distributed  mov¬ 
ing  sources,  and  dynamic  turbulence  fields.  In  this  pa¬ 
per,  we  summarized  recent  efforts  in  developing  a  rigorous 
FDTD  implementation  of  sound  propagation  in  a  moving, 
inhomogeneous  fluid.  As  an  illustration,  calculations  of 
low-frequency  sound  propagation  through  an  atmospheric 
large-eddy  simulation  were  presented.  We  have  also  suc¬ 
ceeded  in  formulating  and  numerically  implementing  new 
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Fig.  11:  Comparison  of  the  exact  impedance  equation 
from  the  relaxation  model  Eq.  (7)  (solid  lines)  and  its  ap¬ 
proximation  Eq.  (8)  (dashed  lines). 
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Fig.  12:  Propagation  above  a  soft  porous  ground  in  the 
presence  of  a  rigid  barrier,  0.250  s  after  source  initiation. 
Calculation  was  performed  with  a  time-domain  boundary 
condition  for  the  ground  interface. 


time-domain  methods  for  handling  sound  interactions  with 
the  ground  surface.  The  methods  are  unique  in  that  they 
are  based  on  full-frequency  equations  for  sound  propaga¬ 
tion  in  porous  materials  including  viscous  and  thermal  dif¬ 
fusion  processes. 

Taken  together,  the  progress  reported  here  lays  the 
foundation  for  highly  realistic  simulations  of  atmospheric 
sound  propagation.  Among  the  applications  are  virtual 
testing  of  tactical  acoustic  sensors  and  their  associated  in¬ 
formation  processing  algorithms.  We  anticipate  that  this 
new  capability,  when  applied  to  sensor  design  and  procure¬ 
ment  for  the  Army’s  Future  Force,  will  substantially  lower 
development  costs  and  time  cycles,  and  ultimately  improve 
sensor  performance  in  complex  battlefield  environments. 
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